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In  the  present  paper  we  describe  an  algorithm  for  the  evaluation  of  Bessel  functions 
Y„(x)  and  Hi^\x)  {j  =  1,2)  of  arbitrary  positive  orders  and  arguments  at  a  constant  CPU 
time.  The  algorithm  employs  Taylor  series,  the  Debye  asymptotic  expansions  and  numerical 
evaluation  of  the  Sommerfeld  integral,  and  is  based  on  the  following  two  observations. 

1)  The  Debye  asymptotic  expansions,  contrary  to  what  appears  to  be  a  popular  belief, 
are  not  expansions  in  inverse  powers  of  (large)  parameter  v  but  turn  out  to  be  uniform 
expansions  in  inverse  powers  of  (large)  parameter  g\  =  {x  —  v)lx^l^  for  i  >  i/  and  (large) 
parameter  g2  =  {y  —  x)lv^l^  for  x  <v. 

2)  For  X  and  v  such  that  both  Taylor  and  Debye  expansions  do  not  provide  a  specified  accu¬ 
racy  Bessel  functions  can  be  computed  at  a  constant  CPU  time  via  (numerical)  evaluation 
of  the  Sommerfeld  integral  along  contours  of  steepest  descents. 

In  addition,  in  Appendix  B  we  obtain  certain  new  estimates  concerning  decay  of  the  func¬ 
tions  Ju{x)  and  — l/y„(a:)  of  fixed  x  and  large  i/,  and  in  Appendix  C  we  show  that  functions 
Ju{x)  of  integer  v  provide  the  solution  for  a  certain  system  of  coupled  harmonic  oscillators. 
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1.  Introduction 


Bessel  functions  of  argument  x  and  order  v  of  the  first  kind  J^(x),  second  kind  V’„(x)  and 
third  kind  ni^\x),  hI^\x)  (Hankel  functions),  play  an  important  role  in  physics,  mathematics 
and  engineering.  Applications  of  Bessel  functions  usually  require  an  algorithm  for  the  rapid 
evaluation  of  these  functions  with  sufficiently  high  accuracy. 

For  arguments  x  ~  1  and  arbitrary  t/  Bessel  functions  can  be  computed  via  their  Taylor 
expansions  (see  Subsection  2.2  below).  If  i  >  1  and  i/  <  x5  these  functions  can  be  evaluated 
by  means  of  the  Hankel  asymptotic  expansion  (see  Subsection  2.3  below).  However,  there  exists 
a  wide  region  of  values  of  x  and  v  where  both  Taylor  and  Hankel  expansions  do  not  provide 
any  reasonable  numerical  approximation. 

Most  of  the  existing  algorithms  for  the  evaluation  of  Bessel  functions  in  this  region  are 
based  on  the  recurrence  relation  (see,  for  example  [1]) 

(1) 

where  f^(x)  denotes  one  of  the  functions  J„(x),  Yi,{x)  or  ny\x)  (j=l,2).  The  asymptotic 
estimate  of  the  complexity  of  these  algorithms  is  of  order  0(x)  for  functions  of  the  first  kind 
and  of  order  0{y)  for  functions  of  the  second  and  third  kind  (see,  for  example,  [1]). 

In  this  paper  we  present  an  algorithm  for  the  evaluation  of  an  individual  Bessel  function  of 
an  arbitrary  nonnegative  order  and  argument  at  a  constant  CPU  time.  The  method  is  based 
on  the  following  two  observations. 

1)  The  Debye  asymptotic  expansions  [3],  proposed  in  1909  and  since  that  time  considered 


expansions  for  large  orders  (see,  for  example  [1],  [4],  [5],  [6]),  are  found  to  have  a  much  wider 
range  of  validity.  Namely,  we  show  that  for  x  >  z/  this  expansion  for  function  H^\x)  is  a 
uniform  asymptotic  expansion  in  inverse  powers  of  (large)  parameter  gi  =  (x  —  i')/x3  (see 
Theorem  3.1  and  Observation  3.1  below).  Moreover,  it  turns  out  that  for  i/  =  0  the  Debye 
asymptotic  expansion  coincides  with  the  Hankel  asymptotic  expansion  (see  Theorem  3.2  below). 
For  X  <  V  the  Debye  expansions  for  functions  J^(x)  and  Yu{x)  axe  proved  to  be  uniform 
asymptotic  expansions  in  inverse  powers  of  (large)  parameter  g2  =  (i'  —  x)/i/^  (see  Theorem 
3.3  and  Observation  3.2  below). 
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2)  For  the  values  of  x  and  v  for  which  both  Taylor  and  Debye  expansions  do  not  provide 
a  specified  accuracy,  Bessel  functions  can  be  computed  (at  a  constant  CPU  time)  by  means  of 
numerical  evaluation  of  the  Sommerfeld  integral  taken  along  Debye  contours  (the  definitions  of 
the  Sommerfeld  integral  and  Debye  contours  are  presented  in  Subsection  2.4  below).  It  is  worth 
noting  that  Debye  contours  were  extensively  investigated  in  connection  with  the  derivations 
of  various  asymptotic  expansions  for  Bessel  functions  (see,  for  example,  [3],  Ch.  8  of  [4], 
[7]).  However,  the  possibility  of  using  them  in  numerical  computations  seems  to  have  been 
overlooked. 

The  plan  of  the  paper  is  as  foUows.  In  Section  2  we  summarize  certain  mathematical 
facts  to  be  used  in  the  rest  of  the  paper.  In  Section  3  we  analyze  the  error  terms  of  the  Debye 
asymptotic  expansions.  Numerical  evaluation  of  the  Sommerfeld  integral  is  discussed  in  Section 
4.  In  Section  5  we  briefly  discuss  the  implementation  of  our  numerical  scheme.  In  Section  6 
we  present  a  formal  description  of  the  algorithm.  In  Appendix  A  we  discuss  round-off  errors. 
In  addition,  in  Appendix  B  we  obtain  asymptotic  solutions  (with  respect  to  u)  of  equations 
Ji,{x)  =  €  and  -l/Yt,{x)  =  e  for  large  fixed  positive  x  <  u  and  sufficiently  small  c  >  0.  Finally, 
in  Appendix  C  we  show  that  functions  Jt,{x)  of  integer  i/  describe  displacements  of  coupled 
harmonic  oscillators  on  a  line. 

2.  Relevant  mathematical  facts 

In  this  section  we  present  a  number  of  well  known  formulae  to  be  used  in  the  rest  of  the  paper. 

2.1  Connections  between  the  three  kinds  of  Bessel  functions 

All  the  formulae  presented  in  this  subsection  can  be  found,  for  example,  in  [1]. 

Functions  Hu\x)  {j  =  1,2)  are  expressed  through  J„{x)  and  Yv{x)  as 


H’i^\x)  =  J^{x)  +  i  Y,ix), 

(2) 

Hl^\x)  =  Mx)  -  i  n(x). 

(3) 

For  nonnegative  x  and  i/,  both  Jv{x)  and  Y^{x)  are  real.  Thus 

Mx)  =  Re  H^ix), 

(4) 
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Yt,(x)  =  Jm  hI^\x), 

(5) 

hI^\x)  =  HI'\x)  -  2ilm  Hl^\x). 

(6) 

2.2  Taylor  expansions 


For  small  arguments,  function  Ju{x)  is  normally  evaluated  via  the  formula 

k 


If  I'  is  not  an  integer,  function  Y^(x)  is  computed  as 
T/  Mx)  cos(7r  i/)  - 

= - — - r - . 

sin(7r  v) 

For  integer  v  the  formula  (8)  can  not  be  used  and  is  replaced  by 

k\  ^ 


(7) 


(8) 


(D'SK)’ 

where  ^(z)  =  ^ln(r(^)).  Formulae  (7)-(9)  can  be  found,  for  example,  in  [1]. 

2.3  The  Hankel  asymptotic  expansion 

The  Hankel  asymptotic  expansion  has  the  form  (see,  for  example,  [1],  Ch.  7  of  [4]  ,  Ch.  7  of 

[5]) 

exp(ix(x,i'))  (”“)  ,  (10) 


where 

x  1 

X(x,  «^)  =  X  -  -  - -XI/, 

and  ^7v+i,a(»'»®)  is  the  error  term. 

Coefficients  bn(i/)  satisfy  the  recurrence  relation 

6o(i/)  =  1, 


(11) 


(12) 
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(13) 


^  ,  ,  (2n  +  1)2  -  4  1/2 

bn+i{v)= - ^  - b„{v),  (n=0,l,---), 

whereas  the  error  term  Oi^+i^h(u,x)  is  bounded  by  (see,  for  example,  Ch.  7  of  [5]) 

|^N+i,ft(t^,j)|  <  2  -  exp  ^  .  (14) 


2.4  The  Sommerfeld  integral 

All  the  formulae,  presented  in  this  subsection  can  be  found,  for  example,  in  Ch.  8  of  [4]. 

For  the  values  of  x  and  v  for  which  both  Taylor  and  Dabye  expansions  do  not  provMe  a 
specified  accuracy  we  computed  function  Hi^\x)  by  means  of  numerical  evaluation  of  the  so 
called  Sommerfeld  integral: 


,  ,  1  foo+t-r 

HI  ’(x)  =  —  I  exp(a;  sinh(«;)  —  uw)  dw. 
irt  J-00 


(15) 


Following  [4]  we  wiU  write 


w  =  u  +  iv,  (16) 

where  both  u  and  v  are  real.  Integration  in  (15)  is  performed  along  an  arbitrary  contour  that 
has  the  following  asymptotes: 

lim  u  =  0,  (17) 

u— *— 00  ' 

Um  V  =  TT.  (18) 

ti^OO  ^  •' 


Observation  2.1 

As  paths  of  integration  in  ( 15)  it  is  natural  to  choose  the  so  called  Debye  contours  on  which  the 
integrand  of  (15)  does  not  oscillate  (see,  for  example,  Ch.  8  of  [4]).  We  note  in  passing  that  the 
Debye  contours  are  a  particular  example  of  contours  of  steepest  descents  that  are  widely  used 
for  the  evaluation  of  the  asymptotic  expansions  of  certain  contour  integrals  (see,  for  example, 
Ch.  8of[4],[5],  [6]). 
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The  Debye  contours  u  =  u{v)  are  curves  on  the  complex  u)-plane  (16)  that  (generally 
speaking)  are  implicitly  defined  via  equations 


For  X  >  v, 


where 


p{u,  v)  =  0. 


p(u,  u)  =  cosh(u) 


cos(P)  = 


sin(/?)  +  {v-  /3)  cos(/?) 

sin(v)  ’ 


It  immediately  follows  from  (21)  that  for  nonnegative  x  and  i/, 


For  X  <  V, 


and 


0</s<|. 


p{u,v)  =  V,  if  u  <  a, 


where 


p(u,  v)  =  cosh(u)  -  cosh(a)  .  .  if  w  >  a, 

sm(T?) 


cosh(a) = 


For  X  =  V, 


and 


p(u,  v)  =  V,  if  u  <  0, 


V 

p(«,  v)  =  cosh(u) - r-r-T?  if  w  >  0. 

sm(u) 


(19) 

(20) 

(21) 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 
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Finally  we  note  that  aU  the  Debye  contours,  associated  with  function  (of  nonnega¬ 

tive  X  and  positive  v)  lie  within  the  strip 

0  <  V  <  X.  (28) 

Graphs  of  Debye  contours  can  be  found,  for  example,  in  Ch.  8  of  [4]. 


2.5  The  Debye  asymptotic  expansions 

In  this  subsection  we  present  formulae  for  the  Debye  asymptotic  expansions  that  can  be  found 
(in  a  slightly  different  form)  in  Ch.  10  of  [5]. 

For  u  <  X, 

Hl'Kx)  =  (lY +  (29) 

Vxy  (i2_j,2)4  U  J 

where 


/  2  2\-  /  ^ 
r]i  =  (x  -  1/  p  -  u  axccos  (  -  1  —  — , 


^N+i,2(i',p)  is  the  error  term,  polynomials  Un{t)  are  defined  in  (35),  (36)  below,  and 

V 


la:2-i/2|5 


For  X  <  u, 


J^{x) 

n(x) 


Zx)2  (l/2  -  i2)4  V  J 


1  +  ^iv+i,i(i'5  0)  (27 


where 


and  ^N+i,i(*')P)  and  9n^\,2{v,p)  are  the  error  terms. 
Polynomials  «„(<)  are  defined  by  the  formulae 


uo(f)  =  1, 


tii(t)  = 


24 


(3t  -  5t^), 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 
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(36) 


«n+l(0  =  ^  (1  ^niT)dr, 

(n  =  0,1,---)- 

The  error  terms  in  the  Debye  expansions  satisfy  the  inequalities 


\^N,2iv,p)\  <  2 exp 

i^o,p(wjv+i) 

I/N+1  ’ 

(37) 

|6n+i,i(«/,p)|  <  2  exp 

J'2Fi,p(«i)' 

\  l^i,p(«Ar+i) 

)  i/^+i  ’ 

(38) 

l^N+i,2(i^,p)l  <  2  exp 

^VoAuiY 

\  1^0,p(WN+l) 

;  i/jv+i  • 

(39) 

In  (37)  -  (39),  symbols  Va,b{f)  and  Vab(/)  denote  the  so  called  total  variations  of  functions 
f(x)  and  f(ix),  respectively  (see,  for  example,  Ch.  1  of  [5]); 

T/  //1A\ 


=  J 

Ja 


3.  Erro*  terms  of  the  Debye  asymptotic  expansions 

In  this  section  we  obtain  estimates  of  the  error  terms  of  the  Debye  asymptotic  expansions  (29), 
(32)  and  (33).  We  start  with  a  more  detailed  analysis  of  the  polynomials  Un(t)  defined  in  (35), 
(36). 

3.1  The  polynomials  Un{t) 

Lemma  3.1 
For  any  n  >  1, 


where 


Unit)  =  t"  Unit), 


^nit)  = 

fc=0 
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The  coefficients  are  defined  by  the  formulae 


aC  =  1, 

aj  =  0  if  A:  <  0  or  fc  >  n,  (n  =  0, 1,  •  •  •), 


,n+l 


=  «!  ( 


n  +  2k 


+ 


J-.  ( 


2  '  8  (2*  +  n  +  1) 

n  +  2A:  —  2  5 


)- 


)■ 


(n  =  0,1,' 


8  (2A:  +  n  +  1) 

A:  =  0, 1,  •  •  •,  n  +  1). 


(44) 

(45) 


(46) 


Proof 

We  will  prove  the  lemma  by  induction.  For  n  =  1  the  formulae  (42)-  (46)  immediately  follow 
from  (35).  Suppose  now  that  the  formulae  (42)-(46)  are  satisfied  for  certain  n  =  m  >  1.  Then 

m+l 

52((m-^2A:)a^*-(m-|-2it-2)a^'_l)•^2^  (47) 

k=o 

and 


/  (1  -  5r^)  Um{T)dT  = 

Jo 

m-H  .. 

y  771 - - - -  5ar-i )  •  (48) 

Now  substituting  (47)  and  (48)  into  (36)  we  observe  that  (42)-(46)  hold  for  n  =  m  -f  1  which 
concludes  the  proof  of  the  lemma.  □ 

The  following  coroUary  is  an  obvious  consequence  of  the  lemma  and  the  formulae  (12),  (13). 
Corollary  3.1 
For  all  n  >  0, 


OO  =  ^n{0), 

where  the  coefficients  6n(t')  are  defined  in  (12),  (13). 


(49) 


8 


Lemma  3.2 
For  any  <  >  >  0, 


«n(»<)  >  Uniitl)  >  0, 

(50) 

«n(iO  ^  «n(*<l)  ;r.  /  «  -_j  .  /  « 

f2n  —  ^2n  ’  if  ^  ^  0  and  ti  ^  0, 

(51) 

dUn{it)  ^  ^ 
dt 

(52) 

V.n{it)  >  ti„(t), 

(53) 

dUn{ii)  ^  dUnit) 

dt  -  dt  ' 

(54) 

dUr^{it)  ^  dUn(t) 

dt  ~dt 

(55) 

Proof 

It  immediately  follows  from  (44)- (46)  that  for  aU  n  >  0  and  k  <  n, 

<  =  {-l)‘aJ,  (56) 

where 

a]?  >  0.  (57) 

Substituting  (56)  and  (57)  into  (43),  we  observe  that  for  any  real  t  all  the  coefficients  of 
the  polynomials  u„(tt)  are  positive  and  therefore  the  inequalities  (50)-(54)  are  satisfied.  The 
inequality  (55)  follows  from  (42),  (53)  and  (54). □ 

Remark  3.1 

While  many  recurrence  relations  occurring  in  mathematical  physics  are  numerically  unstable, 
the  recursion  (46)  is  numerically  stable  since  according  to  (56)  and  (57),  both  terms  in  this 
relation  have  the  same  sign. 

The  following  lemma  is  an  immediate  consequence  of  (42),  (52)  and  (55). 


9 


Lemma  3.3 

For  any  p  >  0  and  n  >  1, 


,p(Wn)  —  f 

Jo 


dt 


dt-p'^-  Un{ip). 


Furthermore, 


dun{t) 

dt 


dt  <  t^O,p(«n). 


3.2  Region  x  >  v 

Theorem  3.1 
For  any 


X  >  u  >Q 


the  error  term  ^Ar+i,2(i',p)  in  the  expansion  (29)  satisfies  the  inequality 


where 


|^N+i,2(i',p)l  <  2 exp 


UN+i{i) 

|(iv+x)’ 

6i 


X  —  V 

9\  =  — i-- 
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Proof 

The  inequality  (60)  and  the  definition  (62)  show  that 


X3 


and  therefore 


Now  combining  (31)  with  (60)  and  (64)  we  observe  that 


1 


(58) 

(59) 

(60) 

(61) 

(62) 

(63) 

(64) 

(65) 
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(66) 


and  substituting  (65)  into  (50)  we  obtain 
Uniip)  <  Un  I 

V 

Combining  the  inequalities  (51),  (63)  and  (66)  we  have 


Un(ip)  <  iin(0  j  ■ 

(67) 

Substitution  of  (31),  (64)  and  (67)  into  (58)  yields 

l'"  (^2  _  i.2)§n 

(68) 

Observing  that 

a.(i)  =  1 

(69) 

and  substituting  (68)  and  (69)  into  (37)  we  immediately  obtain  the  inequality  (61).  □ 

Observation  3.1 

Obviously,  function  ni^\x)  can  be  viewed  not  as  a  function  of  x  and  u,  but  as  a  function 
of  a:  a.  u  the  parameter  gi  defined  in  (62).  Then  the  estimate  (61)  shows  that  for  pi  >  0 
(i.e.  X  >  i/)  the  Debye  asymptotic  expansion  (29)  is  not  an  asymptotic  expansion  in  inverse 
powers  of  (large)  parameter  i/  but  it  turns  out  to  be  a  uniform  (with  respect  to  i)  asymptotic 
expansion  in  inverse  powers  of  (large)  parameter  gi.  Moreover,  as  follows  from  (61),  the  error 
term  0n+i,2{i^,p)  may  be  small  even  if  i/  is  not  large.  The  following  theorem  describes  the 
behavior  of  the  Debye  expansion  (29)  in  the  limit  i/  0. 

Theorem  3.2 

For  any  x  >  0  and  u  =  0  the  Debye  asymptotic  expansion  (29)  and  the  Hankel  asymptotic 
expansion  (10)  are  identical. 

Proof 

From  the  definitions  (11)  and  (30)  we  have 

(^)’  exp(i(x-^)).  (70) 
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Next,  combining  (31)  with  (42)  and  (43)  we  obtain 

Urn  (ip)  =  Hm  ^ ag  =aS(--)”-  (71) 

l/"  "  ^  u-*0  (a.2  _  ^2)§«  \^  -^  )  \  xj 

Now  substituting  (70)  and  (71)  into  (29)  and  taking  into  account  (49)  we  see  that  (for  i/  =  0) 
the  expansions  (10)  and  (29)  are  identical.  □ 

3.3  Region  x  <  v 

Theorem  3.3 
For  any 


z/  >  I  >  0 


(72) 


the  error  terms  ^iv+i,2(*^»p)  of  the  expansion  (32)  satisfies  the  inequality 


\6N+i,2ii',p)\  <  2exp 


\i9i/ 


%+i(0 


92/  92 


§(A/+1)’ 


where 


92 


If  —  X 
i  • 

1/3 


Furthermore,  the  error  term  ^n+i,i(i'»p)  of  the  expansion  (33)  is  hounded  by 


|^A^+i,i(t',p)|  <  2  exp 

Proof 

Substitution  of  (59)  into  (38)  yields 


«Ar+i(i) 


3  ■  P2  /  92 


f  (Af+l)  ’ 


(73) 


(74) 


75) 


(76) 


Now  the  proof  of  the  inequality  (73)  becomes  almost  identical  to  that  of  the  inequality  (61) 
(see  Theorem  3.1)  and  we  omit  it. 

To  prove  (75)  we  observe  that  for  x  <  v, 

V 


P  = 


(v^  - 


r>i- 


(77) 
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Therefore 


Vx,v{Un)- 


\dUn{t) 


dt 


dt  =  Vo,p(ii„). 


Substitution  of  (78)  into  (38)  produces 


(78) 


(79) 


Comparing  (39)  and  (79)  we  see  that  the  proof  of  (75)  is  reduced  to  the  proof  of  (73).  □ 


Observation  3.2 

Parallel  to  Observation  3.1,  we  can  view  function  not  as  a  function  of  i  and  i/  but  as 

a  function  of  i/  and  the  parameter  52  defined  in  (74).  Then  the  estimates  (73)  and  (75)  mean 
that  the  Debye  asymptotic  expansions  (32)  and  (33)  are  not  expansions  in  inverse  powers  of 
(large)  parameter  1/  but  turn  out  to  be  uniform  (with  respect  to  u)  asymptotic  expansions  in 
inverse  powers  of  the  (large)  parameter  52  (compare  with  Observation  3.1). 

Observation  3.3 

In  Appendix  B  it  is  shown  that  for  x  >  1  and  1/  >  x  +  const  •  X3  function  Ju{x)  becomes  small 
whereas  function  lTi/(x)|  becomes  large.  Therefore,  for  sufficiently  large  x  the  most  important 
part  (in  terms  of  applications  of  Bessel  functions)  of  the  region  x  <  u  can  be  estimated  as 

1 

X  <  1/ <  X  +  const  •  x5,  (80) 

Combining  (62),  (74)  and  (80)  we  have 

52  =  IpiI  (1 +  0  (x"3))  .  (81) 

The  estimate  (81)  and  Observation  3.2  show  that  in  the  region  (80)  for  x  >  1  the  Debye 
cisymptotic  expansions  (32)  and  (33)  behave  almost  like  uniform  expansions  in  inverse  powers 
of  (large)  parameter 

V  —  X 

-1-  =  -51-  (82) 
X3 
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4.  Certain  properties  of  the  Sommerfeld  integral 

In  this  section  we  discuss  certain  analytical  properties  of  the  Sommerfeld  integral  (15)  relevant 
to  its  numerical  evaluation. 


4.1  Region  x  >  v 

It  immediately  follows  from  Ch.  8  of  [4]  that  the  explicit  representation  of  the  Sommerfeld 
integral  (15)  on  the  Debye  contours  defined  by  (19)  and  (20)  has  the  form 

hI^\x)  =  ^  exp(i(x  sin(y3)  -  I//?))  exp(i</>i)  du,  (83) 

where 


<l>i  =  <f>i{v,u{v),l3)  =  cos(v)  •  sinh(u)  —  cos(^)  •  u. 


(84) 


and 


du 


(sin(u  —  /?)  —  (v  -  /?)  •  cos(/?)  •  cos(t;)). 


(85) 


dv  sinh(u)  sin^(t;) 

In  (83)-(85)  function  u  =  u{v)  is  evaluated  via  (19),  (20)  and  the  parameter  0  is  defined  in 

(21). 


Theorem  4.1(see  §8.31  of  [4]) 

Function  ^(v,  u(i;),/3),  defined  in  (84),  is  a  nonpositive  decreasing  function  of  \v  —  fi\.  It  has 
the  only  maximum  at  v  =  fi  where 

^(y3,u(/3),/?)  =  0.  (86) 


The  following  corollary  is  an  immediate  consequence  the  theorem. 

Corollary  4.1 
The  equation 

X  =  (87) 

with  €  £  (0, 1),  t;  €  (0,  tt)  and  x  >  0  has  two  and  only  two  solutions  and  /?2  such  that 

fi2.  (88) 
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Theorem  4.2 
For  any  ^  ^  0, 


<Ai  =  -  sin(/3)  (v  -  /3)2  +  0((v  -  /3f). 


(89) 


Proof 

From  (20)  for  any  /3  ^  0  we  have 

u  =  (v-  I3)  +  0{{v  -  /3)^).  (90) 


Substituting  (90)  into  (84)  we  immediately  obtain  (89).  □ 

In  the  rest  of  the  subsection  we  will  estimate  the  domain  on  the  x  —  v  plane  where  the 
integral  (83)  can  be  evaluated  at  a  constant  CPU  time.  We  start  with  the  following  remark. 


Remark  4.1 

For  any  0  <  w  <  1  , 

exp(a:^)  ~  exp  , 

where  function  f{x,i/)  >  0.  This  formula  is  an  immediate  consequence  of  (19)  and  (20). 


(91) 


As  follows  from  (19),  (20),  (84),  (85)  and  Theorems  4.1  and  4.2  the  integrand  of  (83)  is  a 
nonoscillatory  function  of  v.  Moreover,  for  small  \v  — 

(*  +  i)  «p{x*)  = 

exp  (-xsin(|5)(t; -/?)*)  (1  +  6'i(t;  -  ^)),  (92) 

where  ffj(v  —  /?)  =  0(v  —  /3)  is  the  error  term. 

Next,  suppose  that 


and 


X  •sin(/3)  >  1, 


/3> 


Cl 

j  * 

(i  sin(/?))2 


(93) 


(94) 
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where 


Cl  ~  1  (95) 

is  a  (positive)  constant. 

Observation  4.1 

The  condition  (93)  means  that  the  domain  where  the  integrand  of  (83)  is  (numerically)  not 
zero  is  sufficiently  small.  The  condition  (94)  shows  that  the  distance  between  the  maximum  of 
the  integrand  of  (83)  u  =  /?  and  its  singularity  at  v  =  0  (see  Remark  4.1)  is  larger  than  several 
standard  deviations  of  the  Gaussian  in  (92).  Therefore,  if  the  inequalities  (93)  and  (94)  are 
satisfied  then  the  error  term  0\{v  —  /3)  in  (92)  can  be  approximated  (with  high  accuracy)  by  a 
low-degree  polynomial  of  (v  —  /3)  in  the  domain,  where  the  integrand  of  (83)  is  (numerically) 
not  zero.  Thus  in  this  case  the  evaluation  of  the  integral  (83)  by  means  of  the  trapezoidal  rule 
becomes  a  superalgebraicaUy  convergent  procedure.  Moreover,  the  number  of  nodes  of  this 
quadrature  formula  is  (asymptotically  for  large  xsin(/3)  and  Ci)  independent  of  i/  and  x. 

Theorem  4.3 
For  any 

X  >  1  (96) 


and 


V  <  X 


-  Di  X3, 


(97) 


where  Di  ~  1  is  a  (positive)  constant,  the  inequalities  (93)  and  (94)  are  satisfied. 

Proof 

We  will  prove  the  inequality  (93)  first.  From  the  definition  (21)  we  have 


and  thus 


(98) 


X  sin(/?)  =  (x^  —  j  *  . 
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(99) 


Combining  (96),  (97)  and  (99)  we  have 


xsin(/3)  >  (2Z)i)5  •  13  ^1  +  0  (x  3^)  >  l, 


(100) 


which  completes  the  proof  of  (93).  Now  we  turn  to  the  proof  of  the  inequality  (94). 

Owing  to  (22)  we  have  /3  >  sin(/3)  and  therefore  we  can  replace  the  inequality  (94)  by  a 


stronger  one 


sin(/3)  > 


(xsin(/3))2 

Substituting  (98)  into  (101)  we  have 


Using  (97)  we  can  estimate  the  left-hand  side  of  (102)  as 

3 

x5  •  (1  -  >  x5  •  (2jDi  -x"!  -Dl-x'^y  >  i2Di)*. 

Now  it  follows  from  (101)-(103)  that  the  inequality  (94)  is  satisfied  if 


(101) 


(102) 


(103) 


(104) 


where  €  >  0  is  any  (small)  number.  □ 

The  inequalities  (96),  (97)  and  the  estimate  (104)  define  the  domain  on  the  x  -  plain 
where  the  conditions  (93)  and  (94)  are  satisfied  and  therefore  numerical  integration  of  (83) 
by  means  of  the  trapezoidal  formula  can  be  done  at  a  CPU  time  independent  of  x  and  u  (see 
Observation  4.1). 


4.2  Region  x  <  v 

It  follows  from  Ch.  8  of  [4]  that  the  explicit  representation  of  the  Sommerfeld  integral  (15)  on 
the  contours  defined  by  (19),  (23)  is 


1 

I\  =  —  I  exp(x</»2)du, 
TTl  y_oo 


(105) 
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whereas  on  the  contours  (19),  (24)  we  have 

/  exp(a:<^)  (i+^\dv.  (106) 

TTl  Jo  \  dV  J 

Obviously, 

Hi^\x)  =  /i  +  h.  (107) 

In  (105)  and  (106), 

=  <^(u,  o)  =  sinh(u)  -  cosh(a)  ti,  (108) 

4^3  =  ^(w,u(u),q)  =  sinh(u)  cos(t;)  —  cosh(a)  u,  (109) 


and 


du 

dv 


=  cosh(o) 


(sin(t;)  —  V  cos(t;)) 
sin^(v)  sinh(u) 


(110) 


Function  u  —  u{v)  in  (106),  (109)  and  (110)  is  evaluated  via  (19),  (24)  and  the  parameter  a  is 
defined  in  (25). 


Observation  4.2 

Numerical  evaluation  of  the  integral  (105)  can  be  performed  by  means  of  the  Gauss- Legendre 
quadrature  formula  with  the  number  of  nodes  independent  of  v  and  x  because  its  integrand  is 
an  analytic  nonoscillatory  function. 


In  Theorem  4.4  (see  below)  it  is  shown  that  for  any  a  >  0  the  integrand  of  (106)  has 
singularities  in  the  complex  v-plane.  It  turns  out  (see  the  estimate  (111)  below)  that  for  q  <  1 
(  i.e.  v/x  ~  1)  these  singularities  lie  close  to  the  domain  of  integration  in  (106)  which  impedes 
numerical  evaluation  of  this  integral.  Now  we  will  establish  a  domain  on  the  x  —  v  plane 
where  numerical  evaluation  of  (106)  (for  example,  by  means  of  the  Gauss- Legendre  quadrature 
formula)  can  be  performed  at  a  CPU  time  independent  of  x  and  the  estimate  of  this  domain 
is  obtained  in  Theorem  4.7  below. 
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Theorem  4.4 

For  any  a  >  0  the  integrand  o/(106)  has  two  imaginary  complex  conjugate  branch  points  a±. 
Moreover,  for  a  ■<  1, 

a±  =  ±(3)5ai  +  O(a^).  (Ill) 

Proof 

For  imaginary  w  =  is  (s  is  real)  equations  (19),  (24)  become 

cosh(u(is))  =  cosh(Q!)  .  (112) 

sinn(s) 

Therefore,  for  any  a  >  0  there  exists  a  parameter  s  =  sq  such  that 

cosh(u(iso))  =  1  (113) 

and 

0  <  cosh(u(is))  <  1,  |s|  >  So-  (114) 

As  follows  from  (113)  and  (114)  the  points  s±  =  ±so  (or,  equivalently,  the  points  a±  =  ±iso) 
are  the  branch  points  of  the  function  sinh(u)  =  ^cosh^(u)  “  l)*  and  therefore  (see  (106)  and 
(109))  of  the  integrand  of  (106). 

To  prove  (111)  we  must  solve  equation  (113).  Combining  (19),  (24)  and  (113)  we  have 

cosh(a)  =  (115) 

a± 

Expanding  both  parts  of  (115)  in  (small)  parameters  a  and  a±  we  immediately  obtain  estimate 

(111).  □ 

Theorem  4.5  (see  §8.31  of  [4]). 

Function  (f>z  from  (109)  on  the  contours  defined  by  (19),  (24)  is  a  monotonically  decreasing 
function  of  v  with  the  only  maximum  at  v  =  0. 

Theorem  4.6 
For  any  a  71^  0, 

<i>2  =  sinh(a)  —  acosh(a)  -  ^  sinh(a)  +  0{v*). 


19 


(116) 


Proof 

For  small  v  and  o  5^  0  we  have  from  (19),  (24) 

u  =  i  coth(a)t;^  +  0{v*). 

6 

Substituting  (117)  into  (109)  we  immediately  obtain  (116).  □ 


(117) 


As  follows  from  (19),  (24),  (109),  (110)  and  Theorems  4.5  and  4.6  the  integrand  of  (106)  is 
a  nonosciUatory  function  of  v.  Moreover,  for  small  v. 


exp(x<^)  =  exp(— a:  •  (acosh(a)  -  sinh(a))) 
(^1  ^  *)  ^-^a:sinh(Q:)w2^  (1  +  ^2(v)), 


(118) 


where  =  0(v)  is  the  error  term. 

Observation  4.3 

The  local  behavior  of  the  integrands  of  (106)  and  (83)  is  essentially  the  same  (compare  (118) 
and  (92)).  Therefore  the  conditions  under  which  the  integral  (106)  can  be  numerically  evaluated 
at  a  CPU  time,  independent  of  x  and  i/,  are  equivalent  to  the  conditions  (93)  and  (94)  (see 
Observation  4.1).  These  conditions  are 


X  sinh(a)  »  1, 


(119) 


|a±|>2i - 

(isinh(a))2 

where  the  parameters  a±  are  estimated  in  (111)  and 
C2  ~  1 


(120) 


(121) 


is  a  (positive)  constant.  The  condition  (119)  means  that  the  domain  where  the  integrand  of 
(106)  is  (numerically)  not  zero  is  sufficiently  small.  The  condition  (120)  shows  that  the  distance 
between  the  singularities  of  the  integrand  of  (106)  q±  and  the  real  axis  is  larger  than  several 
standard  deviations  of  the  Gaussian  in  (118). 
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Now  the  proof  of  the  Theorem  4.3  can  be  repeated  almost  verbatim  yielding  the  following 
result: 

Theorem  4.7 
For  any 

X  >  1  (122) 

and 

V  >  X  +  D2  x^ ,  (123) 

where 

£>2  ~  1,  (124) 

is  a  (positive)  constant,  the  inequalities  (119)  and  (121)  are  satisfied. 

4-S  Region  x  »  1/ 

In  order  to  construct  an  algorithm  whose  complexity  does  not  depend  on  x  and  1/  in  the  region 

|i  —  i/|  <  const  •  X3 ,  const  ~  1,  (125) 

(t.e.  when  the  conditions  (97)  and  (123)  are  violated)  we  will  consider  numerical  integration 
of  (15)  along  the  Debye  contours  defined  in  (19),  (26)  and  (27).  We  note  in  passing  that  these 
contours  were  extensively  used  for  the  derivation  of  asymptotic  expansions  of  function  hI^\x) 
for  a:  «  J/  (see,  for  example,  Ch.8  of  [4],  [7]). 

Denoting  the  integral  along  the  contour  (19),  (26)  by  Jj  £ind  the  integral  along  the  contour 


(19),  (27)  by  J2  it  can  be  shown  that 

1 

Ji  =  — :  /  exp(i  sinh('u)  —  vu)du, 
irt  y_oo 

(126) 

1  du 

•^2  =  -:  /  ^x.p{x<l>4)  (i  +  -r)dv, 
irt  Jq  av 

(127) 

and 

=  +  (128) 
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where 


</>4  =  04(r,  u{v),  r)  =  sinh(u)  cos(t;)  —  u  +  t  (u  +  i  v), 
r  =  1  - 

X 

du  (sin(t;)  —  v  cos(t;)) 
dv  sin^(t;)  sinh(u) 

In  (127),  (129)  and  (131)  function  u  =  u{v)  is  evaluated  via  (19),  (27). 

Observation  4.4 

The  integral  (126)  is  merely  the  integral  (105)  with  a  =  0  and  therefore  is  can  be  evaluated  at 
a  CPU  time  independent  of  x  and  u  (see  Observation  4.2). 

Observation  4.5 

Obviously,  the  integrand  of  (127)  is  an  analytic  and  (for  sufficiently  small  |r|)  nonoscillatory 
function  of  v.  Therefore  the  integral  (127)  can  be  computed  (for  example  by  means  of  the 
Gauss-Legendre  quadrature  formula)  at  a  CPU  time  independent  of  i/  and  x.  We  will  show, 
however,  that  for  sufficiently  large  x  and  r  >  0  the  integral  (127)  can  not  be  evaluated  without 
an  unavoidable  round-off  error  (see  Observation  4.6  below).  The  inequalities  (141)  and  (142) 
below  estimate  the  range  of  parameters  x  and  i/  where  this  error  is  small  (see  also  Remark  4.2 
below). 

Theorem  4.8 

For  any  0  <  r  C  1  function  Re<f>4  has  the  maximum  at 

W  =  i(3r)5(l-KO(r)).  (132) 

Moreover, 

4'max  =  ^4(^mox?  w(^mo®)?  ^)  ~  g  T2(l  -f  0(t)).  (133) 


(129) 

(130) 

(131) 
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Proof 

From  equations  (19),  (27)  for  small  v  we  have 


Therefore 


and 


«  =  ^  V  +  — +  0(v^). 

35  45-35 

1  2 

sinh(u)  cos(w)  =  -y  v - +  0(t;®), 

35  5-35 


sinh(u)  cos(v)  —  u  = - +  0(i;®). 

9-35 

Substitution  of  (134)  and  (136)  into  (129)  produces 

Re  4>4  =  \  T  V - ^  +  0(t  v^)  +  O(v^). 

32  9-35 

The  formulae  (132)  and  (133)  are  a  consequence  of  (137).  □ 


Theorem  4.9  (see,  for  example,  Ch.8  of  [4]). 
In  the  region  (125)  for  x  >  1, 

|ff(»(x)|=0(x-5). 


(134) 

(135) 

(136) 

(137) 

(138) 


Observation  4.6 

It  follows  from  the  estimate  (133)  that  for  sufficiently  large  x  the  integrand  of  (127)  may  be 
large  whereas,  according  to  Theorem  4.9,  the  integral  itself  is  (asymptotically)  small.  Therefore 
in  this  case  there  exist  cancellations  that  account  for  the  round-off  errors. 

We  do  not  expect  significant  cancellations  if 

X<t>max  <  C'a,  (139) 

where 

C3~l  (140) 
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is  a  (positive)  constant  and  <i>max  is  estimated  in  (133).  The  condition  (139)  means  that  in 
order  to  avoid  cancellations  the  maximum  of  the  integrand  (127)  must  be  of  order  0(1). 

Using  (130)  and  (133)  it  is  easy  to  show  that  for 

x>l  (141) 

the  condition  (139)  is  equivalent  to  the  condition 

u  <  X  -  Dz  x3  +  0(x~3)^  (142) 

where 

03  =  (3  03)5^1,  (143) 

In  fact,  owing  to  (138),  there  must  be  another  condition  of  the  absence  of  cancellations  in 
addition  to  (139).  Namely,  the  domain  in  v  where  the  integrand  of  (127)  is  not  small  must  be 
of  order  0(x~3).  It  can  be  shown,  however,  that  this  condition  holds  if  (142)  is  satisfied. 

We  will  briefly  discuss  the  case  of  r  <  0.  Now  the  integrand  of  (127)  does  not  have  a 
maximum  for  v  ^  0.  However,  in  the  vicinity  of  w  =  0  this  integrand  is  of  order  0(1).  On  the 
other  hand,  function  J„{x)  becomes  small  if  roughly  speaking  (see  Appendix  B) 

j/ >  X  +  const  X3.  (144) 

Therefore  function  J^(x)  (i.e.  the  real  parts  of  Hi^\x))  can  not  be  evaluated  by  means  of  (127) 
without  an  unavoidable  round-off  error  if  i/  satisfies  (144)  with  const  >  1.  It  can  be  shown, 
however,  that  for 

1/  <  X  +  X>4X3  ,  (145) 

where 

D4  ~  1  (146) 

is  a  (positive)  constant,  both  real  and  imaginary  parts  of  function  hI^\x)  can  be  evaluated 
without  a  significant  round-off  error. 
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Remark  4.2 

The  numerical  computation  of  the  integrals  (126)  and  (127)  provide  a  method  for  the  evaluation 
of  function  Hi^\x)  at  a  constant  CPU  time  in  the  region  (125)  i.e.  where  the  algorithms 
discussed  in  Subsections  4.1  and  4.2  fail  (compare  the  inequalities  (97),  (123)  and  (142),  (145), 
as  well  as  the  estimates  (104),  (124)  and  (143),  (146)).  Moreover,  for  i  ~  1  and  i  ~  i/  the 
integral  (127)  can  be  evaluated  without  a  significant  round-off  error  even  if  (142)  is  violated 
because  in  this  case  the  maximum  of  its  integrand  is  of  order  0(1)  (see  Theorem  4.8  and 
inequality  (139)). 

5.  Implementation  of  the  algorithm 

L.  sec.  ion  we  present  certain  details  of  the  implementation  of  the  algorithm  for  the  evalu¬ 
ation  of  function  .Y^^^(i)  via  the  Debye  asymptotic  expansions  (see  Section  3)  or  the  contour 
integration  (see  Section  4).  This  scheme  was  tested  to  provide  double  precision  accuracy  (at 
least  ti.  rteen  digits)  for 

2  <  I  <  100000  (147) 

and 

0  <  1/  <  100000  -I- 16  •  1000005 .  ( 148) 

For 

z  <  2, 

1/  >  0  (149) 

function  can  be  evaluated  by  means  of  Taylor  expansions  (see  Subsection  2.2).  Discus¬ 
sion  of  the  round-off  errors  for  both  Taylor  expansions  and  the  contour  integration  is  presented 

in  Appendix  A. 

5.J  Implementation  of  the  Debye  asymptotic  expansions 

Formulae  (42)-(46)  make  numerical  evaluation  of  the  Debye  asymptotic  expansions  (29),  (32) 
and  (33)  fairly  straightforward.  However,  we  must  estimate  the  values  of  parameters  x  and  v 
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for  which  these  expansions  provide  a  specified  (in  our  case  double  precision)  accuracy.  It  is 
necessary  to  point  out  that  the  computation  of  a  JV-term  Debye  expansion  involves  evaluation  of 
n-order  polynomials  Un{t)  for  n  =  0, 1,  •  •  • ,  iV  (see  Lemma  3.1).  In  other  words,  the  complexity 
of  this  procedure  is  of  order  0{N'^)  operations  and  therefore  for  given  x  and  v  it  is  desirable 
to  choose  N  as  small  as  possible.  The  following  considerations  provide  a  simple  recipe  for  the 
optimal  (for  given  x  and  u)  truncation  of  the  Debye  expansions. 

The  estimates  (61),  (73),  (75)  and  Observation  3.3  show  that  for  any  fixed  x  >  1  and  given 
N  there  exists  a  constant  pyv  >  1  such  that  the  error  terms  of  the  expansions  (29),  (32)  and 
(33)  (truncated  after  N  terms)  become  small  if  x  and  u  satisfy  the  inequality 

jt/ —  x|  >  x^.  (150) 

Our  numerical  experiments  showed  that  for 

X  >  17  (151) 

and  u  satisfying  (150)  the  Debye  expansions  (29),  (32)  and  (33)  provide  double  precision  ac¬ 
curacy.  The  estimates  of  constants  are  presented  in  Table  1;  these  values  were  obtained  by 
means  of  both  analyzing  the  inequalities  (61),  (73),  (75),  and  comparing  the  estimates  provided 
by  the  Debye  expansions  with  that  by  contour  integration. 

The  optimal  (for  given  x  and  i/)  truncation  of  the  Debye  expansions  can  be  done  by  the 
following  procedure.  First  we  compute  the  parameter  g  =  |x  —  t'l/x? ,  then,  using  Table  1,  find 
gN  <  g  closest  to  g  and  retain  N  terms  corresponding  to  this  g^.  For  example,  if  p  =  8  then 
giv  =  7  and  thus  N  =  13.  As  we  see  from  Table  1  the  Debye  expansions  fail  to  provide  double 
precision  accuracy  if  p  <  6.5. 

Remark  5.1 

In  addition  to  the  error  term  ®Ar+i,i(i',p),  estimated  in  (75),  the  expansion  (32)  contains  the 
error  term  ^a^+i,i(i',  0).  It  can  be  shown,  however,  that  if  the  inequalities  (150)  and  (151)  are 
satisfied  then  (with  double  precision  accuracy)  this  term  can  be  neglected. 
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5.2  Numerical  integration  for  x  >  v 

As  follows  from  (92),  for  laxge  xsin(/3)  the  integrand  of  (83)  is  sharply  peaked  at  v  =  j3,  and 
thus  the  interval  in  v  where  this  integrand  is  not  small  may  be  much  narrower  than  the  actual 
interval  of  integration  v  €  (0,  tt).  To  estimate  numerically  meaningful  domain  of  integration 
in  (83)  we  must  find  (unique)  solutions  of  equation  (87)  /3i  and  /?2  (see  Corollary  4.1)  with 
c  approximately  equal  to  the  absolute  value  of  the  (specified)  error  of  the  evaluation  of  (83). 
After  that  numerical  integration  in  (83)  is  restricted  to  the  intervad  v  6  (/3i,/?2)- 
In  accordance  with  Observation  4.1  and  Theorem  4.3  it  was  found  that  for 


and 


with 


a:  >  10 


V  <  X  -  di 


1 


(152) 


(153) 


d\  «  5,  (154) 

the  minimal  number  of  nodes  of  the  trapezoidal  formula  (see  Observation  4.1)  needed  to  e.al- 
uate  (83)  with  double  precision  accuracy  is 

n  =  25  (155) 

independent  of  z  and  v.  If  (152)  or  (153)  are  violated  then  to  obtain  the  same  accuracy  we 
have  to  increase  n.  Changing  the  variable  of  integration  (see,  for  example,  [8]) 

V  =  t\  (156) 

it  is  possible  to  somewhat  improve  the  estimate  (154).  It  was  found  that  after  the  change  of 
variable  (156)  we  can  evaluate  (83)  with  double  precision  accuracy  in  the  region  (153)  where 

di  w  1.5,  (157) 

with  the  number  of  nodes  (of  the  trapezoidal  formula) 

n  =  37  (158) 
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independent  of  x  and  v. 


5.3  Numerical  integration  for  x  <  v 

We  will  first  consider  the  integral  (105)  (see  Observation  4.2).  Function  ^2?  defined  in  (108) 
on  the  interval  — oo  <  u  <  a  has  the  only  maximum  at  u  =  — q.  Moreover,  for  large  x 
this  function  is  not  small  only  for  u  «  —a.  Therefore,  we  must  first  estimate  the  interval 
u  €  (o!i,a!2)  (a  ,  —a  <  02  <  a)  where  the  integrand  of  (105)  is  (with  given  accuracy)  not 
zero  and  restrict  the  numerical  integration  to  this  interval.  We  evaluated  (105)  by  means  of 
the  Gauss-Legendre  quadrature  formula.  It  was  found  that  for 

X  >  1  (159) 

the  minimal  number  of  nodes  required  to  evaluate  this  integral  with  double  precision  accuracy 
is 

n  =  45  (160) 

independent  of  x  and  u. 

Now  we  will  discuss  the  computation  of  the  integral  (106).  As  follows  from  (118)  for 
sufficiently  large  x  sinh(Q)  the  integrand  of  (106)  is  sharply  peaked  at  v  =  0.  Therefore  we 
must  first  estimate  the  interval  v  €  (0,03)  (03  <  jt)  where  the  integrand  of  (106)  is  (with 
given  precision)  not  zero  and  restrict  the  numerical  integration  to  this  interval  (compare  with 
Subsection  5.1). 

The  integral  (106)  was  evaluated  by  means  of  the  Gauss-Legendre  quadrature  formula. 
In  accordance  with  Observation  4.3  and  Theorem  4.7  it  was  found  that  this  integral  can  be 
evaluated  with  double  precision  accuracy  with  the  number  of  nodes 

n  =  45  (161) 

independent  of  x  and  v,  if  x  satisfies  the  inequality  (159)  and 

1/  >  x-\-d2  X3,  (162) 
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where 


da  «  1.5.  (163) 

5.^  Numerical  integration  for  x  ^  v 

The  problem  of  the  numerical  evaluation  of  the  integrals  (126)  and  (127)  is  essentially  the  same 
as  that  of  the  integrals  (105)  and  (106)  (see  Subsection  5.2)  and  here  we  will  briefly  discuss  the 
results  of  our  numerical  experiments. 

The  integrals  (126)  and  (127)  were  evaluated  by  means  of  the  Gauss-Legendre  quadrature 
formula.  In  accordance  with  Observation  4.6  and  the  estimates  (142),  (144)  it  was  found 
that  both  real  and  imaginary  parts  of  these  integrals  can  be  computed  with  double  precision 
accuracy  with  the  number  of  nodes 

n  =  33  (164) 

independent  of  x  and  v,  if 

X  >  10,  (165) 

and 

|z  -  <  ds  -zs,  (166) 

where 

da  «  2.  (167) 

In  addition  (see  Remark  4.2),  this  scheme  provides  double  precision  accuracy  with  n  defined  in 
(164)  for 

1  <  z  <  10  (168) 

and 

x>v.  (169) 
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6.  Description  of  the  algorithm 

In  this  section  we  describe  the  algorithm  for  the  evaluation  of  an  individual  function  H^\x) 
(and  therefore,  according  to  (4)  -  (6),  all  other  Bessel  functions)  of  nonnegative  x  and  v.  The 
final  algorithm  consists  of  two  parts:  algorithm  A  that  combines  contour  integration  and  Taylor 
expansion,  and  algorithm  B  that  is  based  on  the  Debye  asymptotic  expansions. 

6.1  Algorithm  A 

if  X  <  2  then  compute  H^\x)  via  Taylor  expansions  (see  Subsection  2.2). 
end  if 

if  X  >  2  and  v  >  x  +  1.5-X3  then  compute  H^\x)  by  means  of  evaluation  of  the  Sommerfeld 
integral  along  contours  (19),  (23),  (24)  (see  Subsections  4.2  and  5.3). 
end  if 

if  X  >  9  and  v  <  x-1.5-x3  then  compute  n^\x)  by  means  of  evaluation  of  the  Sommerfeld 
integral  along  contours  (19),  (20).  (see  Subsections  4.1  and  5.2). 
end  if 

if  X  <  9  or  \y  -  x\  <  1.5  •  x?  then  compute  n^\x)  by  means  of  evaluation  of  the  Sommerfeld 
integral  along  contours  (19),  (26),  (27)  (see  Subsections  4.3  and  5.4). 
end  if 

6.2  Algorithm  B 

if  X  >  17  and  i/  <  x  —  6.5  •  X3  then  compute  H^\x)  by  means  of  evaluation  of  the  Debye 
expansion  (29)  (see  Subsections  4.1  and  5.1). 
end  if 

if  X  >  17  and  i/  >  x  +  6.5  •  X3  then  compute  Hi^\x)  by  means  of  evaluation  of  the  Debye 
expansions  (32),  (33)  (see  Subsections  4.1  and  5.1). 
end  if 

6.3  The  final  algorithm 

if  X  >  17  and  \x-u\>  6.5  •  xi  then  compute  by  means  of  the  algorithm  B 
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else  compute  by  means  of  the  algorithm  A 

end  if 

6.4  Timing 

A  computer  (FORTRAN)  program  using  the  algorithm  described  in  the  preceding  sections  was 
implemented  and  tested  on  Sun  SPARCstation  1.  This  program  consists  of  approximately  2000 
executable  lines. 

We  compared  the  time  for  the  evaluation  of  function  hI^\x)  by  our  algorithm  with  the  time 
required  to  compute  this  function  by  means  of  recurrence  relation  (1)  (for  increasing  orders) 
in  case  of  integer  1/.  It  was  found  that  in  the  range  of  validity  of  the  Debye  expansions  our 
algorithm  catches  up  with  the  recursion  for  1/  «  10;  in  the  region  where  the  contour  integration 
is  used  the  same  happens  for  1/  «  800.  For  arguments  x  <  17  the  algorithm  is  approximately 
20  times  slower  than  the  recursion. 

b.  Conclusions 

In  the  present  paper  we  have  shown  that  the  Debye  asymptotic  expansions,  contrary  to  what 
appears  to  be  a  popular  belief,  are  not  expansions  in  inverse  powers  of  (large)  parameter  v 
but  turn  out  to  be  uniform  expansions  in  inverse  powers  of  (large)  parameter  =  (i  -  v)/x3 
foT  X  >  u  and  (large)  parameter  g2  =  {v  —  x)/vi  ioi  x  <  v  (see  Theorems  3.1,  3.3  and 
Observations  3.1  and  3.3).  For  x  and  v  such  that  both  Taylor  and  Debye  expansions  do  not 
provide  a  specified  accuracy  we  have  demonstrated  that  function  H^\x)  can  be  computed 
at  a  constant  CPU  time  via  (numerical)  evaluation  of  the  Sommerfeld  integral  along  contours 
of  steepest  descents  (the  Debye  contours).  Obviously,  numerical  integration  along  contours  of 
steepest  descents  can  be  applied  for  the  evaluation  of  other  functions  of  mathematical  physics. 
In  particular,  it  can  be  used  for  the  computation  of  Bessel  functions  of  complex  arguments  and 
orders  (the  classification  of  the  Debye  contours  in  case  of  complex  x  and  v  can  be  found,  for 
example,  in  Ch.  8  of  [4]). 

In  addition,  we  have  obtained  new  estimates  concerning  decay  of  functions  Ju{x)  and 
-l/y„(x)  of  fixed  X  >  0  and  large  positive  v  (see  Appendix  B).  Finally,  we  have  shown  that 
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Bessel  functions  of  the  first  kind  and  integer  orders  provide  a  solution  to  a  system  of  differential 
equations  for  a  chain  of  coupled  harmonic  oscillators  (see  Appendix  C). 

The  author  would  like  to  thank  Professor  Vladimir  Rokhlin  for  useful  discussions  and  for 
his  interest  and  support. 
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Appendix  A 
Round-off  errors 

In  this  Appendix  we  briefly  discuss  round-off  errors  that  appear  (for  certain  values  of  z  and  v) 
when  function  hI^\x)  is  evaluated  via  either  Taylor  expansions  or  contour  integration. 

A.l  Taylor  expansions 

It  is  weU  known  (see,  for  example,  Ch.  3  of  [4])  that  for 

u  =  m  a,  (170) 

where  m  is  an  integer  and  \<t\  ■<  1, 

J^{x)  cos(7r  I/)  a  J_^(x),  (171) 

and  thus  in  this  case  formula  (8)  produces  significant  round-off  error.  However,  for  any  fixed 
z  >  0  function  yi,(z)  is  an  analytic  function  of  i/  and  therefore  it  can  be  evaluated  by  means 
of  interpolation  with  respect  to  the  order. 

Our  experiments  showed  that  that  for  |<t|  >  5  •  10“^  no  significant  error  occurred.  For 
smaller  |cr|  we  used  Chebychev  interpolation  of  function  Fl/(z)  on  the  interval 

1/ €  [m  -  10"^  m -I- 10“^]  (172) 

with  the  number  of  nodes  k=6;  this  number  of  nodes  proves  to  be  sufficient  for  the  evaluation 
of  this  function  with  at  least  thirteen  digits  for  u  from  (172)  and  z  <  2. 

A. 2  Numerical  integration  for  x  >  v 

Our  experiments  showed  that  for  z  >  1  it  is  impossible  to  compute  the  integral  (83)  without 
a  round-off  error  unless  function  defined  in  (84),  is  carefully  evaluated.  This  error  occurs 
because,  as  follows  from  (83),  for  z  >  1  small  errors  of  the  evaluation  of  function  <t>i  produce 
large  errors  of  the  integrand  of  (83). 
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fi{v,l3)  =  {v- 13)  cos{/3)  (175) 

hiv,  (3)  =  sm{v)  -  sm(/3)  (176) 

we  see  that  for  small  |r  —  /3|  each  of  the  functions  fi{v,/3)  and  f2{v,f3)  axe  of  order  0{v  -  (3) 
whereas  the  numerator  of  (174)  is  of  order  0({v  — 

To  avoid  the  round-off  error,  cased  by  the  cancellation  of  the  leading  terms  of  the  Taylor 
expansions  of  functions  (175)  and  (176),  we  can  first  evaluate  the  numerator  of  (174)  using 
its  Taylor  expansion  in  (small  parameter)  (v  —  /?)  and  after  that  compute  sinh(u)  by  means  of 
(173),  and  u  via 

u  =  ln(sinh(«)  -|-  cosh(«)).  (177) 

In  addition,  a  round-off  error  appears  if  we  evaluate  function  <f>i  itself  via  the  formula  (84) 
Rewriting  (84)  in  an  equivalent  form 

(pi  =  (cos(t;)  -  cos(/3))  sinh(ii)  -|-  (sinh(u)  —  u)  cos{P),  (178) 

we  see,  that  in  order  to  avoid  cancellations  (and  therefore  the  lost  of  significant  digits  in  (178)) 
we  can  evaluate  functions  (cos(r)  -  cos(^))  and  (sinh(u)  -  u)  via  their  Taylor  expansions  in 
(small  parameters)  (v  -  /?)  and  u,  consequently. 

Our  experiments  showed  that  the  integral  (83)  can  be  computed  without  significant  round¬ 
off  error  if  we  use  Taylor  expansions  in  (174)  and  (178)  for  |v  —  /?|  <  0.1. 


A. 3  Numerical  integration  for  x  ^  v 

It  was  found  that  the  integral  (126)  can  be  computed  without  round-off  error  for  any  x  and 
V.  However  it  turns  out  that  we  cannot  evaluate  of  the  integral  (127)  without  a  round-off 
error  which  becomes  laxge  for  a:  >  1.  Like  in  case  x  >  v.  the  main  source  of  this  error  is 
the  sensitivity  of  the  computation  of  the  integrand  of  (127)  (for  large  x)  to  small  errors  of  the 
evaluation  of  function  <f>4  defined  in  (129).  To  analyze  this  effect  we  observe  that  as  follows 
from  (134)  and  (135)  for  small  v  functions  u  and  sinh(u)  •  cos(u)  are  of  order  0{v)  whereas 
their  difference  (136)  is  of  order  O(u^).  Therefore  the  round-off  error  of  the  evaluation  of  (136) 
(and  thus  of  function  04  from  (129))  appears  (for  sufficiently  small  v)  due  to  cancellation  of 
the  leading  terms  of  the  Taylor  expansions  of  (134)  and  (135). 

It  follows  from  our  numerical  experiments  that  if  the  left-hand  sides  of  (134)-(136)  are 
evaluated  via  their  Tailor  expansions  for  v  <  1  then  the  integral  (127)  can  be  computed 
without  significant  round-off  error.  These  expansions  are: 

u  =  0.57735026918962576  v  -f  0.02566001 1963983367  -I- 

0.0014662863979419067  A  0.000097752426529460445  + 

0.74525058224720925  •  10~®  u®  -|-  0.61544207267743328  •  10"®  + 

0.5290.118464628039 •  10"’'  +  (179) 

sinh(u)  =  0.57735026918962576  v  +  0.057735026918962576  -I- 

0.0062775386411887880  v®  -I-  0.000655246734  8028954  u’’  -b 
0.000066970892258993254  v®  -f  0.67971226373232793  •  10"®  -|- 

0.68878126140038453  •  10"®  -I-  •  •  • ,  (180) 

u  -  sinh(«)  cos(u)  =  0.25660011963983367  u®  -I- 

0.00097752426529460445  u’’  -b 
0.000072409204836637368  u®  + 

0.74478039260541289  •  10"®  u”  -f 
0.74130822294291681  •  10"®  +  ■  ■  ■ .  (181) 
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Appendix  B 

Decay  of  functions  J„(a;)  and  — l/y^(i)  for  the  large  orders  and  fixed  arguments 

In  this  appendix  we  discuss  the  behavior  of  functions  J^(i)  and  Y^{x)  for  fixed  arguments  and 
large  orders.  Theorems  B1  and  B2  and  formulae  (219)  and  (222)  contain  the  principal  results 
of  this  appendix. 

Bl.  Statement  of  the  problem 

It  is  well  known  that  for  any  fixed  i  >  0  and  f  oo  function  J^(x)  decays  rapidly  (see,  for 
example,  Ch.  10  of  [5])  and  has  an  asymptotic  behavior  defined  by  the  formula 

However,  this  approximation  is  numerically  meaningful  only  for  i/  >  x^/4  (see,  for  example, 
Ch.  10  of  [5]).  On  the  other  hand,  often  it  is  necessary  to  have  an  accurate  estimate  of  an 
order  i/j  >  x  >  0  such  that 

J^(x)  <  €,  (183) 

for  all  1/  >  i/j.  In  (183)  i  is  fixed  and  €  >  0  is  supposed  to  be  sufficiently  small.  This  problem 
arises,  for  example,  when  one  implements  Miller’s  algorithm  (  see,  for  example,  [1],  [2]),  or 
sums  a  Neumann  series 

00 

=  X]  '^M+n(a:).  (184) 

n=0 

It  is  well  known  that  the  behavior  of  function  7„(x)  of  fixed  positive  arguments  and  large 
orders  is  close  to  that  of  function  —  l/y„(i).  For  example,  for  any  fixed  i  >  0  and  i/  — >  oo  this 
function  decays  rapidly  (see,  for  example,  [1])  and  has  the  following  asymptotic  representation 

In  this  appendix  we  prove  that  functions  J„(x)  and  -l/yj,(x)  of  any  fixed  x  >  0  are 
monotonically  decreasing  functions  of  i/  >  x  (see  Theorems  Bl  and  B2  below).  In  addition, 
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in  case  of  large  fixed  x  >  0  and  sufficiently  small  c  >  0  we  present  approximate  solutions  of 
equations 


Deforming  the  contour  of  integration  in  (192)  into  the  contour  (19),  (24)  with  v  €  (— Tr,^)  we 
obtain 


dJ,{x) 
du 


I  r ,  ^  f  du\ 

=  ~2^  J  y  ^  sinh(ti)  cos(t;)  -  i/u)dv,  (193) 

where  the  derivative  dujdv  is  defined  in  (110).  On  this  contour, 


u  >  0, 

du{v)  _  du(—v) 
dv  dv 


(194) 

(195) 


Combination  of  (193),  (194)  and  (195)  yields 


dJ„{x) 

dv 


1 

ir 


u  +  V 


du 

dv 


^  exp(isinh(u)  cos(t;) 


vu)dv. 


Finally,  for  v  €  (0,7r)  we  have  from  (110) 


du 


The  conclusion  of  the  theorem  follows  from  (194),  (196)  and  (197).  □ 


(196) 


(197) 


Observation  Bl 

It  follows  from  formula  (188)  and  Theorem  Bl  that  for  any  v  >  x  >  0  (x  is  fixed)  equation 
(186)  has  at  most  one  real  solution.  Moreover,  if  such  a  solution  exists,  then  the  inequality 
(183)  is  satisfied  for  all  r'  >  vj. 


B3.  Certain  properties  of  function  Vvix)  for  v  >  x. 


Turning  to  the  discussion  of  the  behavior  of  function  Yiy{x)  for  j/  >  i  >  0  we  will  first  prove 
the  following 


Lemma  Bl 

For  any  v  >  x  >  0, 


y.(x)  <  0. 


(198) 
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Proof 

Combining  formulae  (5),  (105)  (106)  and  (107)  we  have 

1  /•£»  \  du 

Y^(x)  =  —  /  exp{x4>2)du - /  ex-p{x4>3) —dv,  (199) 

TT  QQ  X  Jo  uV 

where  the  parameter  a  and  functions  <j>2  and  <^3  are  defined  in  (25),  (108)  and  (109),  respectively, 
and  function  u  =  ti(t;)  is  evaluated  via  (19),  (23).  Now  the  conclusion  of  the  lemma  follows 
from  (197)  and  (199).  □ 

We  will  now  prove  the  analogue  of  Theorem  B1  for  function  Y^(x). 


Theorem  B2 
For  any  1/  >  x  >  0, 


dY^jx) 

dv 


<  0. 


(200) 


Proof 


We  stcJt  with  Nicholson’s  formula  (see,  for  example,  Ch.  9  of  [5]): 

N„{x)  =  -^  /  Ko{2x  sinh(t))  •  cosh.{2ut)dt  =  J^{x)  +  Y^{x), 

in 


where 


(201) 


Kq{z)  =  f  exp(-2:cosh(t))df  (202) 

Jo 

is  Macdonald’s  function  of  zero  order.  It  immediately  follows  from  (201)  and  (202)  that  for 
any  x  >  0, 


dNu{x) 

du 


>  0. 


(203) 


Next,  formula  (201)  yields 


^n(x)  _  1  fdN,{x)  _  ,_,dJ.{x)\ 

du  2Y^(x)  V  5*/  du) 


(204) 


The  conclusion  of  the  theorem  is  a  consequence  of  (188),  (190),  (198),  (203)  and  (204).  □ 


The  following  observation  is  closely  related  to  Observation  B1  above. 
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Observation  B2 

Lemma  Bl  and  Theorem  B2  show  that  for  any  fixed  positive  x  <  1/  equation  (187)  has  at  most 
one  real  solution.  Moreover,  if  such  a  solution  exists, 


1 


(205) 


for  aU  i'  >  Vy. 

In  the  rest  of  the  appendix  we  derive  approximate  (asymptotic)  solutions  of  equations  (186) 
and  (187). 

B4.  Asymptotic  solution  of  equation  (186). 

We  will  first  prove  the  following  simple 
Lemma  B2 
For  any  \z\  <  1, 

00  2n-l  1 

+  +  (206) 

n=l  ^ 


Proof 


Expanding  left-hand  side  of  (206)  into  Taylor  series  we  have 

00  ,n  00  -2n-l  1  00  _2n 

n=l  n=l  n=l 

which  concludes  the  proof  of  the  lemma.  □ 


(207) 


An  approximation  to  vj  from  (186),  hereafter  denoted  by  Vj,  will  be  sought  as  a  solution  of  the 
equation 


exp(-^)  ^  ^ 

(27r)2  {i/j  —  x^)* 


(208) 


where  the  function  on  the  left-hand  side  of  (208)  is  the  leading  term  of  the  Debye  asymptotic 
expansion  (32)  with  the  change  of  notation  u  Vj. 
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Introducing  the  notation 


Vj  = 


_ 


(209) 


and  substituting  (209)  into  (34)  we  obtain 

ri2  =  x-yj-  [ln(2/j)  +  In  (1  +  (1  -  y~^)^)  -  (1  -  yj^)^  ■ 


(210) 


n=l 


Noticing  that  for  2  =  (1  —  ^)5  equation  (206)  becomes 

In  (1  +  (1  -  yj^)^)  =  (1  -  yj^)^  +  2^^^^  ■  yj^)^  -  Hyj),  (211) 

and  substituting  (211)  into  (210)  we  obtain 

2n  +  1 


n2  =  x  yj  (1  -  2/7^)2  Yi 


(212) 


n=l 


Now,  substituting  (209)  and  (212)  into  (208),  we  have 

x-yj-(l-  yf)^  Y  2n  +  l  4  ^  ((27r®)2  e)  .  (213) 

nsl 

Finally,  introducing  a  new  (unknown)  function  qj  which  is  a  positive  solution  of  the  equation 

(214) 

and  substituting  (214)  into  (213)  we  have 


yj  =  x  3  gj  +  1, 


3 

T  ^  i  In  gj  +  i  ln(l  +  g?)  = 


20 


—  ln(2*  TTS  13  e).  (215) 


We  seek  the  asymptotic  solution  of  (215)  under  the  condition 

X  >  -ln(2«  ttJ  xi  i)  >  1.  (216) 

Taking  into  account  (216),  equation  (215)  immediately  yields  the  leading  term  of  the  asymptotic 
expansion  of  qy. 


Qi  ~ 


(217) 
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where 


Sj  =  33  2  ^  ln(2^  ar5  . 


(218) 


Corrections  to  (217)  can  be  found  by  standard  methods  (see,  for  example,  Ch.  1  of  [5]).  After 
some  algebra  we  obtain 

^  -  1  =  (<,  i-i)'  (1  -  g  ln(«,)  S-^  +  i  {ij  x-iy  + 

O  (ln(«j)  6f)  +  0  (x-i  6-^))  .  (219) 


B5.  Asymptotic  solution  of  equation  (187). 


In  this  subsection  we  briefly  discuss  the  derivation  of  the  asymptotic  solution  of  equation 
(187);  this  approximate  solution  will  be  denoted  by  i7y.  Like  in  case  of  equation  (186)  we  seek 
an  approximation  to  i/y  as  a  solution  (for  fixed  positive  x  <  i')  of  equation 


(220) 


where  the  function  on  the  left-hand  side  of  (220)  is  the  modulus  of  the  inverse  of  the  leading  term 
of  the  Debye  asymptotic  expansion  (33)  with  the  chajage  of  notation  v  Vy.  The  technique 
employed  for  solving  equation  (220)  is  the  same  as  that  for  equation  (208)  and  we  omit  the 
computational  details.  It  can  be  shown  that  the  asymptotic  solution  of  (220)  derived  under 
the  condition 


X  >  —  ln(2 


1 

4X2 


X  3  e)  >  1, 


(221) 


has  the  form 

^-1  =  («y  ln(5j,)^;3  +  ^(^^aj-i)" + 

o  (in(«,)  ef)  -h  o  (x-l  e;^)) ,  (222) 

where 

6y  =  35  2“2  ln(2«  7r“2  x”5  c)^®  .  (223) 
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B6.  Comparison  of  exact  and  asymptotic  solutions  of  equations  (186)  and  (187). 

Tables  2  and  3  contain  approximations  j7j  and  t7y,  obtained  via  (219)  and  (222),  as  well  as  exact 
numerical  solutions  of  equations  (186)  and  (187)  for  several  values  of  x  and  e  (when  solving 
equations  (186)  and  (187)  their  left-hand  sides  were  computed  via  contour  integration  for  a:  >  2 
and  Taylor  expansion  for  x  <  2).  It  is  interesting  to  note,  that  the  formulae  (219)  and  (222) 
provide  reasonable  approximations  to  i/j  and  i/y  even  for  x  «  1,  i.e.  when  the  conditions  (216), 
(221)  are  violated. 

Remark  B1 

It  is  easy  to  see  from  (219)  and  (222),  that  up  to  logarithmic  (in  x)  corrections  the  parameters 
i7j  and  i7y  can  be  estimated  as 

a  X  -f  Cj(e)  X3,  (224) 

z7j,  «  x-i-cy(e)  X3,  (225) 

where  Cj(€)  «  Cy(€)  >  0  are  independent  of  x.  In  other  words,  the  approximations  (224),  (225) 
provide  a  (rough)  estimate  of  a  domain  on  the  x  —  i/  plane  where  functions  J„(x)  and  -  1/Fi,(x) 
of  any  fixed  x  >  1  are  small  for  all  >  Vj  a  Vy. 

Appendix  C 

Bessel  functions  and  a  chain  of  harmonic  oscillators 

In  this  appendix  we  show  that  functions  J^(x)  of  integer  v  describe  displacements  of  coupled 
harmonic  oscillators  on  a  line. 

Differentiating  the  formula  (see,  for  example,  [1]) 

^  =  5(/.-i(*)  -  /.+■(*)).  (226) 

where  f„{x)  has  the  same  meaning  as  in  formula  (1)  we  obtain 

-  2  Ux)  +  /,«(!)).  (227) 
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Equation  (227)  can  be  compared  with  the  differential  equations  for  a  system  of  equal  point 
masses  on  a  line  which  interact  with  their  nearest  neighbors  via  the  elastic  force.  The  displace¬ 
ment  from  the  equilibrium  Un(t)  (as  a  function  of  time  t)  of  an  n-th  such  oscillator  satisfies 
equation  (see,  for  example,  [9]) 

M  =  G  (u„_i(t)  -  2  «„(t)  -b  u„+x(f)),  (228) 

where  M  is  the  mass  and  G  is  the  elastic  constant. 

The  analogy  between  (227)  and  (228)  becomes  especially  transparent  if  we  rewrite  (227) 
for  Bessel  functions  of  the  first  kind  of  integer  order  v  =  n.  It  follows  from  (227)  that  these 
functions  of  even  orders  satisfy  the  system  of  differential  equations 

=  4(«^2n-2(l)  -  2  J2n(x)  -h  (u  =  0,  ±1,  ±2,  •  •  •);  (229) 

with  the  initial  conditions  (see,  for  example,  [1]) 

Jo(0)  =  1,  </2n(0)  =  0,  (n  =  ±1,  ±2,  •  •  •); 

=  0,  (n  =  0,±l,±2,---).  (230) 

x=0 

For  odd  orders  we  have  the  same  system  of  equations 

- ^^2'^  ^  ~  4(«^2n-i(®)  - 2  J2n+i(a;)  +  </2n-(-3(a^))5  («  =  0,  ±1,  ±2,  •  •  •);  (231) 

but  the  initial  conditions  in  this  case  are  different  (see,  for  example,  [1]): 


dhnjx) 

dx 


^2n+i(0)  =  0,  (n  =  0,±l,±2,---); 


dJi{x) 
dx  ,=o 
d</2n+l(®) 
dx 


dJ-i(x) 

dx 


(n=  1,±2,---). 


(232) 


Comparing  (228)  with  (229)  and  (231)  we  see  that  function  J2n{x)  (or,  equivalently,  J2n-i-i(a;)) 
can  be  viewed  as  a  displacement  of  an  n-th  oscillator  at  a  ’time’  a:  in  a  chain  (228)  with 


parameters 


M  =  1, 

G  =  i.  (233) 
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It  is  interesting  to  observe,  that  as  follows  from  (229)  and  (231)  the  chains  of  even  and  odd 
Bessel  function  do  not  interact. 

This  analogy  enables  us  to  give  a  mechanical  interpretation  of  certain  properties  of  Bessel 
functions. 


1.  A  zero  of  a  Bessel  function  can  be  interpreted  as  a  ’time’  at  which  a  corresponding 
oscUlator  passes  the  equilibrium.  Therefore,  the  well  known  result,  that  function  Jn{x)  has 
infinitely  many  zeros  as  a:  ^  oo  reflects  the  obvious  physical  property  that  any  oscillator  (in  a 
chain  with  zero  friction)  passes  the  equilibrium  infinitely  many  times  (as  time  goes  to  infinity). 

2.  The  identity  (see,  for  example,  [1]) 

f;  J^(x)  =  J^ix)  +  2  f;  J^(x)  =  1  (234) 

n=— oo  n=l 

means  that  the  oscillators  in  the  chains  (229)  and  (231)  under  the  initial  conditions  (230)  and 
(232)  oscillate  in  such  a  way  that  the  sum  of  squares  of  the  displacements  in  both  of  them  does 
not  depend  on  ’time’  i. 

3.  In  terms  of  the  mechanical  model  the  approximation  (219)  (or  its  simplified  version 
(224))  estimates  the  range  of  propagation  of  the  initial  perturbation  (230)  (or  (232))  at  the 
’time’  X. 


4.  It  is  easy  to  show  that  the  well  known  relation  (see,  for  example,  [1]) 
'^2n(®)  =  M^)  +  2  ^2n(x)  =  1 


(235) 


n=— oo  n=l 

is  a  consequence  of  the  conservation  of  momentum  in  the  even  chain  (229).  To  prove  this  we 
observe  that  according  to  the  initial  conditions  (230)  the  total  initial  momentum  of  the  even 
chain  is  zero.  Because  this  system  is  isolated  we  can  write 


dJ2n(x)  _  _ 
dx 


(236) 


for  any  ’time’  x.  Integrating  (236)  and  observing  that  due  to  (230)  the  constant  of  integration 
is  equal  to  unity  we  immediately  obtain  (235). 
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5.  An  interesting  formula  can  be  obtained  from  the  law  of  conservation  of  energy  in  the  odd 
ciiain  (231).  Using  (233)  we  find  the  kinetic  K  and  the  potential  Ft  energy  of  the  odd  chain 
(231)  (see,  for  example,  [9]): 


n  =  i  E 

ns=— oo 


(237) 

(238) 


Combining  (226)  and  (238)  we  obtain 


(239) 


As  follows  from  (232)  and  (233)  the  total  initial  energy  of  the  odd  chain  is  equal  to  1  /4  and, 
the  system  being  isolated,  it  remains  the  same  at  any  ’time’  x.  Therefore  from  (237)  and  (239) 
we  have 


(240) 


The  formula  (240)  means  that  the  sum  of  kinetic  energies  of  the  chains  (229)  zmd  (231)  is 
independent  of  ’time’  x. 


Remark  Cl 

As  follows  from  the  preceding  analysis  Bessel  functions  are  a  (rare)  example  of  a  discrete 
dynamic  interacting  system  where  the  coordinate  of  any  particle  can  be  computed  at  a  CPU 
time  independent  of  the  physical  time. 
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Table  1:  Numerical  estimates  of  the  parameter  gs- 


N 

5 

9 

13 

17 

9N 

23 

10 

7 

6.5 

Table  2;  Comparison  of  the  numerical  solution  of  (186)  with  the  approximation  (219). 


€ 

ESI 

I  =  10 

I  =  50 

X  =  100 

X  =  1000 

X  =  10000 

X  =  100000 

t'o 

6.25 

8.37 

20.10 

66.38 

1040.8 

10082.1 

100163.9 

D 

isi 

6.63 

8.61 

20.22 

66.44 

■Ml 

1040.8 

10082.1 

100164.0 

1^0 

10.28 

13.15 

78.79 

135.8 

10156.0 

100326.7 

11.20 

13.74 

27.83 

78.88 

135.9 

10156.0 

100326.7 

t'o 

21.20 

39.76 

98.33 

10267.8 

100569.4 

D 

IB 

■rIBI 

22.85 

40.33 

98.53 

10267.8 

100569.4 

t'o 

23.37 

28.34 

50.25 

114.8 

MBSM 

10359.7 

100768.1 

D 

IB 

28.14 

31.47 

51.32 

115.2 

■fawB 

BB 

10359.7 

lj0768.1 

Table  3;  Comparison  of  the  numerical  solution  of  (187)  with  the  approximation  (222). 


c 

x  =  2 

X  =  10 

H 

II 

o 

X  =  100 

X  =  1000 

X  =  10000 

X  =  100000 

a 

m 

Bl 

9.89 

9.88 

23.02 

22.96 

128.2 

128.1 

fciiiM 

10139.0 

10139.0 

100309.4 

100309.4 

l7y 

m 

142.5 

142.4 

BtwB 

10202.1 

10202.1 

100443.3 

100443.3 

IS 

B 

18.31 

20.59 

22.57 

23.96 

lii 

165.7 

165.8 

10305.2 

10305.2 

100663.6 

100663.6 

a 

B 

m 

29.68 

32.56 

52.43 

53.26 

118.7 

118.9 

185.6 

185.6 

1181.9 

1181.9 

10392.9 

10392.9 

100851.4 

100851.4 
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